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In recent experiments on ultracold matter, molecules have been produced from ultracold atoms by 
photoassociation, Feshbach resonances, and three-body recombination. The created molecules are 
translationally cold, but vibrationally highly excited. This will eventually lead them to be lost from 
the trap due to collisions. We propose shaped laser pulses to transfer these highly excited molecules 
to their ground vibrational level. Optimal control theory is employed to find the light field that will 
carry out this task with minimum intensity. We present results for the sodium dimer. The final 
, target can be reached to within 99% if the initial guess field is physically motivated. We find that 

■ the optimal fields contain the transition frequencies required by a good Franck-Condon pumping 

' scheme. The analysis is able to identify the ranges of intensity and pulse duration which are able 

, to achieve this task before other competing process take place. Such a scheme could produce stable 

^ . ultracold molecular samples or even stable molecular Bose-Einstein condensates. 

D 
0^ 



I. INTRODUCTION 



\^ [ The creation of cold molecules from atomic Bose-Einstein condensates (BEC) as well as 

\^ . from ultracold thermal gases ^ J^] has advanced remarkably over the last two years. In both cases molecules are 

■ formed due to interaction of atoms with an external field. The latter can be the electric field of a laser leading to 
C^l , photoassociation or a magnetic field tuned to drive the atoms across a Feshbach resonance Tuning close 
CD ■ to a Feshbach resonance can furthermore be used to obtain a very large scattering length. This enhances three- 

' body recombination and molecules are thereby formed 0. In most experimental schemes translationally cold, but 

■ vibrationally very highly excited molecules are produced. These highly excited molecules are not stable with respect 
' to collisions, and dimers consisting of bosonic atoms in particular are very rapidly. In the case of sodium this happens 

within milliseconds ^ . Creation is therefore only a first step toward novel experiments using ultracold molecules , 
and their stabilization is the obvious next task. 

The ultimate goal is to create u = 0, ./ = molecules. The task of transferring the highly excited molecules to 
w = 0, J = is far from being trivial: a long-range molecule and a, v — molecule are very different (Cf. Fig. ^ where 
the wavefunctions for the last bound and the ground vibrational level of the ground state of Na2 are drawn) . This 
!• ■ task has to be solved under three major constraints: 1.) The goal has to be achieved in a time short compared to the 
.5^ ' collisional lifetime. The rate for collisional decay is not precisely known at present, but lifetimes on the order of or 
shorter than ms should be expected @. Calculations for collisional decay of molecules in highly excited vibrational 
$H levels have been performed in the non-reactive ^'"^He -I- H2(w) case |l4j . showing an increase of the relaxation rate by 
I nearly three orders of magnitude when the initial state is going from v=l (rate ~ 5x 10~^^ cm^ s~^) to v=10. In the 
reactive Na -I- Na2(w) case, quantum calculations for the v = 1 ^ rate yield ^ 5 x 10^^*'cm'^ s""'^ Therefore, 
at an atomic density of 10^^ cm~'^, a vibrational quenching relaxation time well below 1 ms is indeed to be expected. 
2.) If laser pulses are to be used, spontaneous emission has to be avoided. The radiative lifetime of the excited state 
is on the order of 10 ns. 3.) Due to vibrational energy "pooling" or vibration-to-vibration ladder climbing [l^ the 
intermediate range of binding energies needs to be avoided at any cost, i.e. the molecules need to be immediately 
transferred to the lowest levels. 

In this paper, we suggest to employ optimal control to transfer the highly excited molecules to the rovibrational 
ground state. Optimal control has been intensely studied both theoretically and experimentally in many areas of 
physical chemistry. Optimal control theory (OCT) |l7j offers the prospect of driving an atomic or molecular system 
to an arbitrary, desired state due to the interaction with an external field. Experimentally, control is achieved via 
feedback learning loops, see e.g. for a recent review. The difference between the approach we suggest and control 
experiments as they have been performed over the last decade lies in the different target: In the latter the goal consists 
usually in varying the ratio between different dissociation channels, i.e. many final states are available and the target 
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is rather broad, while in the present context a single ro-vibrational level is to be achieved. However, this does not 
pose a problem in principle - it might make control harder to achieve, but it does not render it impossible. We 
suggest optimal control to transfer the highly excited molecules to the rovibrational ground state because it provides 
an efficient method operating on a very short timescale. 

Alternatively, two-photon Raman transitions have been suggested Q to transfer the highly excited molecules to 
low-lying vibrational states. However, if really the lowest vibrational levels shall be populated, a Raman scheme with 
continuous wave (CW) laser pulses is not conceivable due to unfavorable Franck-Condon overlaps. The framework of 
Stimulated Raman Adiabatic Passage (STIRAP) is not any more promising since a highly oscillatory wave function 
with 80% or more of its probability density in the asymptotic region can not be transferred adiabatically into a very 
compact wave function, typically of Gaussian shape (Cf. Fig.^). This has been overlooked in previous studies [l9ll20| 
in which the vibrational degree of freedom was not explicitly taken into account. 

Optimal control theory has been applied before in the context of the formation of ultracold molecules, specifically 
to calculate optimal conversion of an atomic into a molecular BEG In this work, a laser pulse was employed 

together with a magnetic field which induced a Feshbach resonance. The authors worked in the framework of the Gross- 
Pitaevski equation which is a non-linear Schrodinger equation while employing a linear variant of OGT. A nonlinear 
version of OGT as required for the application to the Gross-Pitaevski equation has been worked out recently ,2^]. The 
vibrational structure of the molecules was simply taken into account by assigning a few vibrational levels in Ref. |2l| , 
but the coordinate-dependence of the vibrational wave functions, in particular the qualitatively different character 
of initial and final state wave functions, was neglected. Finally, the initial state in Ref. corresponds to two free, 
colliding atoms. If one was to include the vibrational structure in the model of Ref. the resulting Hilbert space 
would be infinite-dimensional, and due to the theorem of controllability [23Ll2^ one would not be guaranteed that an 
optimal solution exists at all. 

For conceptual clarity, we will therefore separate the creation of molecules, however weakly bound, from the process 
of their stabilization. This approach corresponds furthermore nicely to the current status of experiments on cold 
molecules which create weakly bound dimers via magnetic Feshbach resonances or three-body recombination |^ |^ |^ 
IE 1^ Qi IS ' The stabilization of these molecules in the experiment still remains an open problem to which no easy 
answers exist. In this paper, we will hence start from extremely weakly bound molecules and employ optimal control 
theory to obtain short, shaped laser pulses which drive the system from a specified, highly excited vibrational level 
to the ground state. 

We will apply the optimal control algorithm to the formation of stable sodium dimers. From molecular spec- 
troscopy experiments with GW lasers |2^, we know that a pathway connecting the rovibrational ground state with 
the last bound levels exists. Furthermore, sodium has been one of the first systems to be studied by femtosecond 
spectroscopy |2^ and control experiments |23| . Its properties under ultracold conditions are equally well studied, see 
e.g. 5, "2^. The scheme we suggest is similar in spirit to that of Ref. 0| where nanosecond pulses were used. In order 
to compare to current experimental control techniques, however, we employ pulses of femtosecond and picosecond 
timescale. These pulses are optimized by the control algorithm while they would serve as input for a feedback learning 
loop in a prospective experiment, see e. g. This is in contrast to the intuitively chosen pulses of Ref. "2^. A 

large number of theoretical studies on femtosecond pulse control of sodium has been reported, mostly in the context 
of photodissociation and ionization, see e.g. |30l IsiL Is^ . Our approach has been motivated by the experiment of 
Ref. '25II where the last bound levels of the Na2 ground state were probed. This was achieved by exciting v = Na2 
molecules from a molecular beam with GW lasers to an excited electronic state and subsequent spontaneous emission. 
A scheme with multiple excitation-deexcitation steps was necessary due to Franck-Gondon overlaps. In principle, one 
could invert this scheme to create v — from highly excited molecules. However, the yield would be low and the 
required times cannot compete with the time scales of collisional loss and spontaneous emission. We therefore look 
for optimal ultrashort pulses which realize this inverse scheme. 

We formulate the problem as one of optimizing state-to-state transitions. This is in contrast to a density matrix 
formulation which describes transitions from an ensemble of states to another one H^l^l- As an extreme case, the 
latter includes true cooling as a transition from an ensemble of states to a single state. Note, that true cooling requires 
coupling to a dissipative environment which accepts the entropy (35.. .36] . However, the optimization of state-to-state 
transitions should be sufficient for Bose-condensed samples since one starts with a coherent sample. Even in the case 
of thermal samples state-to-state transitions can serve as an important first step in which one state is selected to be 
transferred while all others are ionized. We will employ OGT in the Krotov variant restricting the change in pulse 
energy at each iteration H^iH^H^. This guarantees monotonous and smooth convergence toward the objective. 

We will work in the framework of the linear Schrodinger equation which describes the internuclear dynamics of 
two atoms, i.e. in case of a Bose-Einstein condensate (BEG) we neglect the condensate dynamics. This approach is 
justified by the time scales present in standard optical and/or magnetic traps. While the internuclear dynamics and 
pulse shaping occur on the time scale of femtoseconds to picoseconds, the condensate dynamics for conventional traps 
is characterized by microseconds. The condensate will have to adjust to the new internal state, but this is going to 
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happen on a much slower time scale than the one of the pulse [Sg • We can hence safely ignore the influence of the 
condensate on the dynamics of the state-to-state transition and assume that the adjustment of the condensate takes 
place once the pulse is over. 

The paper is organized as follows: The model for the sodium dimer as well as the proposed scheme to form stable 
ultracold molecules are described in Section ^1 Section UTTI briefly reviews the optimal control algorithm with details 
given in the Appendix. The results are presented in Section IIVI Section concludes with a discussion of the 
experimental feasibility of our proposed scheme as well as its implications to vibrational cooling. 



II. THE NAa SYSTEM 



We have chosen Na2 because this system has been intensely studied in ultracold experiments, in control experiments, 
and by traditional spectroscopy, and a large amount of experimental data is available. In particular, highly accurate 
potential energy curves have been obtained from molecular spectroscopy |39l l40l | . Experiments |25l | have furthermore 
shown that a route from the last bound levels to v = exists. In this scheme, ground state molecules from a molecular 
beam with v = 0, J = are excited by a CW laser (A = 610 nm) to the A-'^E+ excited state {v' — 15). Those molecules 
which decay to the v = 29 level of the ground state, are excited by a second CW laser (A — 540 nm) to excited state 
levels with v' — 100 — 140. A third CW laser (A = 595 nm) probes the transition between these excited state levels 
and the last bound levels {v = 61 — 65) of the ground state. 



A. Proposed scheme for the formation of ultracold, stable molecules 

We envisage the following two-step scheme for the production of ultracold molecules (Fig.^. First, loosely bound 
molecules are created by tuning an external magnetic fleld to sweep an atomic sample across a Feshbach resonance 0, 
Hi ^ IS or by enhanced three-body recombination in the vicinity of a Feshbach resonance . Another possibility is 
given by applying a weak off-resonant CW photoassociation laser. The last level is so extremely loosely bound that 
even a very weak field perturbs it leading to a coupling with the continuum states Second, a shaped laser pulse 
is applied to transfer the highly excited molecules to v = 0. This second step is of our concern in the present study. 
Following the experiment j25j, we propose to employ fields which induce transitions between the Na2 ground state 
(X'^^j ) and the A^Ej excited state. 

One could imagine a one-step scheme where the optimal pulse is used to directly transfer a continuum wave function 
to the ro-vibrational ground state. However, in the case of continuum states the theorem of controllability (23i. 
is not applicable and one is not guaranteed that an optimal solution exists. Furthermore, the results rnay depend on 
how the continuum is modelled. Since experimentally, step (1) seems to become a standard technique j^l^^l^, we 
prefer to restrict the task of the optimal pulse to the transfer of one bound level to another bound level and to attain 
a mathematically well-defined model. 



B. Model and methods 



We are interested in obtaining an estimate of the feasibility of optimal control experiments with ultracold molecules. 
We therefore restrict our approach to a qualitative model with two channels which correspond to the two electronic 
states of Ref. and we neglect hyperfine interaction, spin-orbit coupling and rotational excitation. These effects 
should be included to make quantitative predictions. 

Taking hyperfine interaction into account would amount to studying a model with more than two channels. At 
the binding energies at which cold molecules are created in the current experiments, a single channel description for 
the ground state is not valid |39| . i.e. the wave function contains singlet as well as triplet components. The shaped 
laser pulse does not alter the coherence of this superposition. A multi-channel treatment would therefore result in a 
pulse with parts corresponding to the excitation pathway of singlet, and other parts corresponding to the excitation 
pathway of triplet molecules. However, in this study we want to emphasize the importance of the internal structure 
of the molecules, i.e. the coordinate-dependence of the atom-atom interaction potentials, for the creation of stable 
molecules. We therefore restrict the model to two channels, one for the ground and one for the electronically excited 
state, and we perform calculations for singlet states. This amounts to neglecting the triplet component of the wave 
function. If the triplet component of the wave function is larger than the singlet component, one would have to 
repeat the present calculations using the corresponding triplet potentials. The qualitative principle, however, remains 
unaltered. 



FIG. 1: Scheme for the production of stable ultracold molecules: Vibrationally highly excited molecules are created due to the 
interaction with an external field (1) and axe transferred to the ground state by a shaped laser pulse (2). The wave functions 
of the last bound and the v = levels are drawn as dashed lines, while the excited state vibrational levels which have good 
Franck-Condon overlap with the last level and with the ground level, respectively, are indicated by dotted lines. Obviously, no 
direct transition exists. 

The argument concerning rotational excitation is along the same lines as in the case of hyperfine interaction. 
Allowing for rotational excitation would amount to treating additional channels with the centrifugal barrier for each 
J added to the atom-atom interaction. The additional channels increase the size of the search space, but at the same 
time allow for more pathways, i.e. for more flexibility. It is thus not clear a priori whether the search for optimal 
pulses will be slowed down or sped up. Moreover, if the initial state contains some rotational excitation, one might 
speculate that this facilitates control since such a wave function would be more bound than one with J = 0. Of 
course, such a claim can only be proven by a detailed investigation which is beyond the scope of the present work. 

We solve the radial Schrodinger equation for the vibrational motion of the sodium dimer, 



with R the distance between the two nuclei. The vibrational wave function consists of two components, ip = 



corresponding to the two channels. The Hamiltonian describing two electronic states and the nuclear degree of 
freedom, R, is given by 



ih-ip{R,t) = (\ip{R,t), 



(1) 





(2) 



where T = denotes the kinetic energy operator, the dipole operator and e{t) the electric field. The ground 

state potential \tg describes the state of X^T,g symmetry which is correlated to the (Ss+Ss) asymptote of sodium, 
and the excited state potential Vg describes the state of A^E„ symmetry correlated to the (3s-|-3p) asymptote. The 



5 



ground state potential has been obtained by an analytical fit to high-resolution spectroscopic data '39| , while an RKR 
potential ^4^] has been matched to an asymptotic expansion j43| for the excited state potential. We assume the dipole 
operator to be independent of R, ft = const. While this approximation needs to be improved, it might well serve 
as a first step since the two concerned electronic states do not show avoided crossings. The electric field is taken 
to be real which means that a fixed (but not specified) polarization is assumed. The initial guess field is given by 
e{t) = eo S{t) cos(wt) where S{t) is a Gaussian envelope function or a sequence of Gaussians and u! is the central 
frequency of the pulse. 

The wave function (p{R, t) is represented on a grid employing the mapped Fourier grid method |43j | which allows for 
accurately representing even the last bound levels of the ground state potential with a comparatively small number 
of grid points {N = 1024). The mapped Fourier grid method introduces, however, unphysical states which can lead 
to spurious effects in the dynamics. For the sodium system, these "ghost" states can be eliminated by using a basis 
of sine functions (instead of plane waves) Q and by choosing the density of points sufficiently high. The mapped 
grid is calculated from the envelope potential of both Vg and Vg, i.e. the same grid is used on both channels. 

The time-dependent Schrodinger equation, Eq. is solved formally by 

ip{R,t)^U{t,0)ipiR,0). (3) 

We employ a Chebychev expansion [ZHl li^ of the time evolution operator, (j{t), which is numerically exact for a 
time-independent Hamiltonian. In our case of an explicitly time-dependent Hamiltonian, Eq. (O is second order in 
the time step. Note, that we cannot use a less expensive propagation method such as the split propagator since the 
grid mapping introduces multiplications between functions in Fourier space and functions in real space, i.e. terms of 
the kind f{R)^. 

The initial state as well as the target state are taken to be eigenstates of the ground state potential. The eigenstates 
are computed simply by calculating a matrix representation of the Hamiltonian in the sine basis and subsequent 
diagonahzation |4J|. 



III. THE OPTIMAL CONTROL ALGORITHM 



We formulate the optimal control for state-to-state transitions, i.e. we want to find the field which drives a specified 
initial state \(pi) to a specified target state \(pf) at the final time t = T. The objective we want to reach can therefore 
be defined as the overlap between the initial state propagated from time t ^ to t — T with field e{t) and the target 
state, 

F=|((p,|U+(r,0;£)|^/)p. (4) 

U"''(T', 0;e) denotes the evolution operator which completely specifies the system dynamics. In the calculations pre- 
sented below, \(pi) will be a highly excited vibrational level of the Na2 ground state, while \(pf) will be chosen to be 
the vibrational ground state of the Na2 ground state. The length of the optimization time interval, T, is an external 
parameter of the calculation. If one compares to experiment, it is related to the bandwidth of the pulse and the 
number of pixels of the pulse shaper. Taking T to be twice the longest vibrational period which occurs in the problem 
is usually a good guess. In the case of long-range molecules, very different time scales are involved which leads to a 
problem of feasibility. We present a work-around in Sec. IIV Gl A field is optimal if it almost completely transfers the 
initial state \(pi) to the target \ipf), i.e. if F attains a value close to one. 

The objective is a functional of the field e{t). However, it depends explicitly only on the final time T. To use 
information from the dynamics at intermediate times, i.e. within the time interval (0,T), we define a new functional 
J, 

J = -F+ [ gis,^)dt, (5) 
Jo 

where the integral term denotes additional constraints over the system evolution. The optimal field will be found 
by minimization of J. The additional constraints provide the connection between the dynamics and reaching the 
objective, they therefore define how the control is accomplished 22] . This becomes particularly apparent in the 
Krotov variant of optimal control theory (see j22l | for a review) which we will employ in the following. The Krotov 
method allows us to derive an iterative procedure which maximizes the original functional F (minimizes ~F) at the 
final time T while modifying the field at intermediate times and guaranteeing monotonic convergence. Details can be 
found in Appendix 1X1 
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The choice of g{e, ip) is not arbitrary, it needs to comply with the requirement of monotonic convergence pl\ (see 
also Appendix fX|l . In this study, we use 

g{e,^)=g{e) = — [e{t)-e{t)]\ (6) 

where e{t) is the field of the previous iteration. S{t) denotes the constraint to smoothly switch the field on and off 
with shape function S{t). We take it to be Gaussian or S{t) = sin^(7rf/T). Eq.(|SJl has the physical interpretation of 
restricting the change in pulse energy at each iteration. It allows for a smooth convergence toward the objective since 
the change in field vanishes when the optimal field is approached ■ The optimization strategy can be controlled by 
the parameter a: A small value results in a small weight of the additional constraint and allows for large modifications 
in the field, while a large value of a represents a conservative search strategy allowing only small modifications of the 
field at each iteration. The parameter a is not directly related to the number of iterations required to find a solution, 
and the best search strategy is usually given by running the optimal control algorithm first with a small and then 
with a large value of a 47]: A small a will excite many pathways, then the solution can be refined with large a. 
With this choice of g{e, tp) the change in field, Ae = e^™ — e°''^, at every time t, < t < T, is given by 



Ae{t) = 

a 



(^.|U"(T,0;e°i'^) \pf) (^/| U"(i, T; e^i'i) A 0(i, 0; e^-) |^.) 
, ■' ^ ' . ■' 

forward backward forward 



(7) 



(see Appendix ^ for the derivation) . The change in field is determined by the " expectation value" of the dipole 
operator ft, i.e. by matching at time t the initial wave function which has been propagated forward in time from t = 
to t with the target wave function which has been propagated backward in time from T to t. Note, that this is not 
an ordinary quantum mechanical expectation value due to the different fields, 6°'*^ of the previous iteration and e""^^ 
of the current iteration. We consider two channels/electronic states for the Na2 system. The wave function therefore 

contains two components, ^^^^ ' ^^"^ dipole operator, , induces electronic transitions. Note that while 

for the initial and target state pe = 0, the propagated states have a non-zero (/Se-component due to the field. 

Eq. Q implies that information from the dynamics in both the past and the future is used to calculate the new field. 
The iterative algorithm is thus defined: We first pick a guess field e°^'^, propagate with U(T, 0; e°^'^) (forward) and 
evaluate the objective. We then propagate the target state \ipf) from time t = T to time t = (backward), storing 
the propagated wave function for all intermediate times. In a third step we obtain the new field e"-™ at each time t by 
evaluating Eq. Q and we use e"™ to propagate \(pi) (forward). These steps have to be repeated until the objective 
has reached the desired value close to one. 

Note that since Ae = e"™ - Eq. O is imphcit in the new field e"''*. A simple, but sufficient remedy to avoid 
solving this implicit equation is found by employing two different grids in the time discretization 37J. The time grid 
for the wave functions has Nt + 1 points and is defined from t = to t = T, while the grid where the field is evaluated 
has Nt points and ranges from t = At/2 to t — T — At/2 with At — T/Nt the spacing for both grids. The new 
field e"™ at the first interleaved grid point t = At/2 is then calculated from the wave function at t = 0, i.e. from 
0(0, 0; e"™) \ipi) = \ipi), while the wave function is propagated from t = to t = At using the field e"'="(Ai/2). This 
process is repeated for all other time grid points. 



IV. RESULTS AND DISCUSSION 



In the proposed scheme, the ultracold molecules are formed in one of the last bound levels, close to the dissociation 
limit. We have investigated states from the whole range of the vibrational spectrum {v = 10, v = 40, v — 62) in order 
to gain a deeper understanding of the physical mechanism behind the control. Recall that the last bound level of the 
ground state of the sodium dimer has quantum number d = 65. 

Optimal pulses transferring all initial states to the target state {v — 0) were found if the guess pulse was chosen 
in a physically sensible way. This means that the pulse spectrum should coincide with the transition frequencies 
between vibrational levels in the ground and excited state which have a good Franck-Condon overlap. The transition 
frequencies with non-negligible Franck-Condon overlap range from about 6000 cm^^ to about 13000 cm^^. This 
corresponds to a transform limited pulse of 5 fs full width half maximum (FWHM). However, not the whole spectral 
range of 7000 cm~^ needs to be covered. The Franck-Condon overlaps of the initial and the target state with the 
vibrational levels of the excited state indicate one or two spectral regions which should be important in the transfer. 
It was sufficient to accordingly choose one or two central frequencies of the guess pulse. Convergence to reach 99% of 
the objective could then be achieved with a relatively small number of iterations. 




FIG. 2: The Franck-Condon overlaps of the intial state {v = 10, 40, 62, dotted lines) and of the target state (i; = 0, dashed 
lines) with all excited state vibrational levels («') are shown vs. frequency. Also shown are the spectra of obtained optimal 
fields (solid black lines) and of the corresponding guess fields (solid grey lines). 



A. Simple search strategies 

In addition to the spectral range covered by the pulse (cf . Fig. ^ , the rate of convergence depends on the two 
control knobs of the optimal control algorithm, the intensity of the pulse or the integrated pulse energy and the length 
T of the optimization time interval. If we choose a comparatively high intensity and long time T {T > T* = 2tt / to = 
2TTh/\Ey — Ey-i\ with Ey the binding energy of level u), we allow the algorithm a lot of freedom. Convergence to 
99% of the objective is then reached very fast (A'it < 30). This corresponds to a few hours of CPU time for v — 10 
and u = 40 and to a few days for w = 62 on a Linux PC (the increase for w = 62 is caused by the necessity to store 
the backward propagated wave functions on disc rather than in memory for long time intervals). 

The resulting pulses are shown in Figs. |21 and |21 A sequence of pulses with 100 fs FWHM was chosen for u = 10 and 
V — 62, while a single pulse of 355 fs FHWM was taken as guess for v = 40. The spectra of the pulses (solid lines in 
Fig. I^J are compared to the Franck-Condon overlaps of the initial (dotted lines) and target (dashed lines) state with 
all vibrational levels of the excited state potential. For the target state w = 0, the Franck-Condon overlaps with all 
excited state vibrational levels are characterized by a simple broad peak between 8700 cm~^ and 11000 cm"-'^. The 
corresponding vibrational levels of the excited state have quantum numbers w' = 1 to w' = 20. 

For initial state v = 10, non-zero Franck-Condon overlaps can be found between 7150 cm~^ and 12500 cm~^ 
corresponding to levels = to v' = 54. In this case, the control task is comparatively simple: There is a number of 
excited state vibrational levels which have good Franck-Condon overlap with both the initial and the target state, i.e. 
a direct pathway exists. The simplest solution is therefore given by a pulse which addresses the transition frequencies 
between these "common" levels [v' = 1 to u' = 20) and initial and target state. This solution is rapidly found by 
the optimal control algorithm: The transition frequencies between these "common" levels [v' — liov' — 20) and the 
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FIG. 3: The optimal fields (black solid lines) corresponding to Fig. |5| Also plotted are the respective guess fields (dashed 
lightgrey lines) and the shape function S{t) (solid darkgrey lines). The differences between optimal and guess fields for v — 62 
can hardly be seen on the time scale of 16 ps, therefore the inset shows the enlargement of a small interval. 



target state are already contained in the guess pulse (grey line in Fig. [2), and the transition frequencies between the 
"common" levels and the initial state around 7800 cm~^ can clearly be seen in the spectrum of the optimal pulse. 
Since the guess pulse is rather intense, other pathways are excited as well. If an initial guess which is even closer to the 
direct pathway solution, i.e. a pulse with two central frequencies is employed, optimal solutions can be found with less 
intensity and larger FHWM, i.e. smaller spectral bandwidth (not shown). Note, that even in this simple case where 
the solution basically can be guessed, the knowledge of only the transition frequencies is not enough. Such a guess 
pulse with two central frequencies transfers less than 1% from the initial to the target state while the time-frequency 
correlations of the optimal pulse allow for more than 99% transfer. 

The situation becomes more involved for higher excited initial states. No direct pathways, i.e. no excited state 
vibrational levels which simultaneously have good Franck-Condon overlap with initial and target state exist. Moreover, 
non-zero Franck-Condon overlaps with the excited state vibrational levels are usually found in different spectral regions 
for the initial and the target state. The guess pulses are therefore chosen to contain two central frequencies addressing 
the relevant transition frequencies. If sufficient intensity is allowed, a solution can rapidly be found. The intensity is 
needed to find other pathways which are not contained in the guess pulse, and the resulting spectrum is rather broad 
(cf. Fig. [3 middle panel). Note furthermore the different time scales in Fig. O where the length of the time interval T 
has been chosen to be at least twice the vibrational period of the initial state. Since the level w = 62 is very close to 
the dissociation limit, its binding energy is very small and its vibrational period is very long. The required T increases 
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up to 240 ps for v = 65. Such long times obviously pose a problem, this will be addressed in Sec. II V CI 

B. Restricting the pulse intensity 

Since optimal control is based on constructive and destructive interference between different quantum pathways, 
the intensity of the pulse must be sufficiently high to excite several such pathways. Furthermore the guess pulse must 
contain sufficient intensity in the relevant modes for the algorithm to find a solution. If the overlap of the initial state 
propagated with the guess field and the target state is very small in Eq. {Tj), the changes in the field are also very 
small. The algorithm then needs a very large number of iterations, or it may not converge at all. Therefore there 
exists a minimum intensity necessary to find a solution. Note that the minimum intensity we have found theoretically 
does not necessarily have to be the same in a corresponding experiment. This is due to the differents ways in which 
a solution is obtained theoretically and experimentally. The relation between the two still needs to be clarified which 
is, however, beyond the scope of the present work. 

The integrated pulse energy is given by 

£p = eocA / \E{t)\^dt (8) 

with E{t) the optimal field, A the area which is covered by the laser {A = irr^ with r = 300/im), c the speed of 
light and eo the electric constant. The integration can be done numerically over the optimal field or analytically over 
the shape function, the latter usually being sufficient to obtain an estimate. Table lists the minimum pulse energy 



V (initial state) 


10 


40 


62 


^puisc 


60 fiJ 


1.5 mJ 


3.9 mJ 



TABLE I; Minimum intensities needed to find a convergent solution transferring the initial state {v = 10, 40, 62) to the ground 
state {v = 0) 

necessary to find an optimal solution within Na < 100 iterations for v = 10, v = 42, and v — 62. The increasing 
difficulty of the control task with increasing vibrational excitation is reflected in the increase in required pulse energy. 

Since part of the intensity is necessary to find the relevant frequencies, one can speculate that it should be possible 
to further reduce the intensity by dividing an optimal field which contains already the relevant frequencies by some 
factor and to restart the optimal control algorithm. This is indeed possible. New optimal solutions within a reasonable 
number of iterations were found after dividing the optimal field of a previous calculation by a factor between two and 
five (reducing the pulse energy by a factor between four and 25). Using an analytical guess with the corresponding 
intensity did not lead to a convergent solution. An example is shown for w = 62 in Fig.^ The spectrum of the weaker 
field becomes broader, in particular new peaks appear at about 6600 cm^^ and 7800 cm^^. The dynamics induced 
by the two optimal fields of Fig. 01 has been analysed by calculating the time-dependent population of vibrational 
levels, \{ipv/v> |(/3g/e(t))p. The number of levels which attain at some time t a population of more than 10% and 5%, 
respectively, are given in Table ^1 The weaker pulse populates a smaller number of vibrational levels. It is interesting 





Eo = 0.01 a.u. 

fpulse = 16 mJ 


Eo = 0.005 a.u. 
£^puisc = 4 mJ 


II rground 
-"'PopXJ.lO 


3 


6 


Ty-ground 
Pop>0.05 


23 


13 


TV rcxcitcd 
-"'PopXJ.lO 


7 


7 


71 rcxcitcd 
^'PopXJ.OS 


24 


16 



TABLE II: Number of vibrational levels of ground and excited state whose population exceeds 10% (A'pop>o.io) and 5% 
(Afpop>o.o5) at some time t {0 <t <T). Note that the 7 excited state levels with more than 10% population are not the same 
for 16 mJ and for 4 mJ. 

to note that the levels which get populated are not the same for the weaker and the stronger pulse. Hence, a different 
solution is found after the reduction of intensity. This is confirmed by the objective F being F = 7.9 • 10~^ after the 
division of the optimal field by a factor of two, but becoming F — 0.97 after 25 more iterations of the algorithm. The 
fact that the field is not simply scaled but a new solution needs to be found, explains why the intensity of the field 
can only be reduced by a factor between two and five, but not by one or two orders of magnitude. Again, while some 
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FIG. 4: The fields (a-b) and spectra (c) of optimal fields for v = 62. The spectrum of the stronger pulse (a) is plotted in black 
while the spectrum of weaker field (b), slightly shifted for visibility, is plotted in grey. The inset shows an enlargement of the 
high-frequency peak. 

information is already contained in the field, a certain minimum intensity is needed to excite other pathways which 
constitute the new solution. 



C. Restricting the optimization time T 



The optimization time T is a parameter of the calculations, and as explained in Sec. IIV Al a good guess is given 
by taking T to be twice the time corresponding to the smallest frequency which occurs in the problem. Since we are 
interested in the last bound levels, this time becomes very large (on the order of nanoseconds for Na2). Such a time 
scale is too long to efficiently compete with collisions and spontaneous emission. 

Simply choosing a smaller value of T does not solve the problem. In this case, the optimal control algorithm did not 
find any solution for many different guess pulses. A comparatively simple remedy to the problem was found by using 
the result of an optimization with a large value of T, restricting it in time, and employing this modified field as initial 
guess pulse for a second run of the optimization algorithm. The restriction in time was done by Fourier-transforming 
the field, deleting points in frequency space and Fourier-transforming it back to the time domain. Keeping only every 
10th value of the field in frequency space reduces the time interval by a factor of 10. The shape function for the 
second optimization run was chosen to be S(t) = sin(7rt/T) where T denotes the reduced time interval. 

Such a reduction of the optimization time has been tested for up to a factor of 20. Remarkably, it did not matter 
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FIG. 5: Convergence of the optimal control algorithm after the total time T (corresponding to the vibrational period of the 
initial state) has been reduced by a factor of 8. 



whether the points in frequency space were deleted symmetrically around w = or not. There is no rigorous upper 
limit to the factor, by which T can be reduced. However, the more T is reduced, the slower becomes the convergence. 
There is furthermore no upper limit to the objective which can be reached with reduced T. This is shown in Fig. [51 
where the objective is plotted vs. the number of iterations of the control algorithm. After reducing T by a factor of 
8, a field which transfers more than 90% of the initial state to w = is found after iVit ^ 1000 iterations. Note, that 
each iteration requires a full time propagtion. Nevertheless, the calculations do not become prohibitively expensive 
due to the shorter propagation time (T/8). The inset of Fig. O shows g[e)dt vs. the number of iterations. The 
integral corresponds to the value by which the objective is increased at each iteration. While this value decreases 
algebraically after about 60 iterations, it is still sufficient to increase the objective from about 20% after A^it = 60 to 
more than 90% after = 1000. The larger number of required iterations is not only due to the restriction in time, 
but also due to a reduction in integrated pulse energy and hence intensity. 

An example of an optimal field and its spectrum after A'it = 1000 are shown in Fig. El The field shows a sequence 
of very short sub-pulses which is characteristic for optimal fields |47|. The spectrum is rather broad, covering about 
9000 cm~^. It should be noted, however, that already the spectrum after the first optimization run with large T was 
rather broad. At this point it is very difficult to decide whether this broadness is physically necessary or or whether 
it is an arteficial by-product of the mathematical algorithm. Unfortunately, it is not possible in a straightforward 
way to include a restriction in bandwidth as condition in the algorithm without ruining the property of monotonous 
convergence (see App.^for a more detailed discussion). 

More insight is gained by analyzing which vibrational levels the optimal field populates during the course of time. 
Fig.Oshows the total population of the two channels (solid black lines) with the optimal field (in grey) which has been 
scaled for comparison. The sub-pulse structure of the field corresponds to population cycling between ground and 
electronically excited state. For the ground state, the population of initial {v = 62) and target {v = 0) vibrational level 
(dashed lines) are also plotted. Finally, the dashed grey lines in Fig. [3 show the sum of populations of all vibrational 
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FIG. 6: The optimal field and its spectrum. The optimization time interval has been reduced by a factor of 8. 



levels, i.e. J2v For the electronically excited state, the solid black and dashed grey lines overlap, i.e. 

only bound levels get excited. This is not surprising: The highest vibrational levels which get populated (Cf. Fig. (HJ 
are those which have good Franck-Condon overlap with the initial state. These levels are still comparatively strongly 
bound {v' = 101 . . . 108 out of about 220 bound levels). For the ground state, the solid black and dashed grey lines 
differ remarkably, corresponding to a non-negligible population of the continuum. However, this does not need to 
worry us. The continuum gets populated at intermediate times for very a short period when the Na2 system is in a 
coherent superposition with the optimal field. The corresponding population is therefore not lost but brought back to 
bound levels at later times and eventually transferred to ?; = 0. Fig. |H1 shows the population of those vibrational levels 
vs. time which acquire more than 5% of the overall population at some time < t < T. Ground state vibrational 
levels are displayed in the upper panel, while excited state vibrational levels are shown in the middle and lower panel. 
The overall number of vibrational levels which get populated is 7 for the ground state {v — 62 and v — have been 
plotted in Fig. [7| for better visibility) and 16 for the electronically excited state. For the ground state (Fig. |Hl(a)-(b)), 
the two neighbouring levels of the initial state get populated at early times. At intermediate times, the population 
is spread over many different levels (and the continuum) while toward the end of the time interval, the population is 
concentrated in three levels {v = 21,22,23). None of this levels is populated for more than 300 fs. This is the time 
which needs to be compared to the time scale of vibrational energy pooling. For the excited state (Fig.|H|(c)-(f)), the 
vibrational levels which get populated can be assigned to three different groups: At early times a number of levels 
{v' = 101 . . . 108) with good Franck-Condon overlap with the initial state get populated. While the population cycles 
back and forth between ground and excited state, these levels are overall populated for a comparatively long time. 
After about 1 ps, intermediate levels {v' = 18,41, 51 . . .53) are populated. Toward the end of the time interval the 
population is concentrated in levels {v' = 7, 8, 9) which have a very good Franck-Condon overlap with the target state, 
V = 0. Also none of the excited state vibrational levels is populated for more than 300-400 fs which should be short 
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FIG. 7: Population analysis of the dynamics with the optimal field for the electronically excited state (a) and the ground state 
(b). The total population |(p/e|^(t))p (solid black lino) is compared to the sum of populations of all feottnii vibrational levels 

(dotted grey line) with the field plotted in the background (lightgrey, scaled for comparison). For the excited 
state, the black and the grey line arc almost identical indicating that the continuum is not populated, while the curves differ 
for the ground state indicating significant population of the continuum at intermediate times. Panel (b) shows furthermore the 
population of initial (v = 62) and target {v = 0) vibrational level (dashed lines). 



enough to compete with vibrational energy pooling. 

Note, that which specific vibrational levels are populated differs for different optimal pulses which are the result of 
initial guess pulses with different frequencies and intensities and of reduction of T by different factors. The overall 
scheme of population, in particular for the excited state with the three groups at early, intermediate and long times, 
is found in all cases. This is not surprising because it simply reflects the fact that a multi-step scheme is needed to 
transfer a very weakly bound state to the vibrational ground state. 



V. CONCLUSIONS 



We have shown that within a two-state model of Na2 optimal fields which transfer highly excited vibrational states 
to the vibrational ground state to more than 99% can be found. We have performed calculations for moderately, 
highly and and extremely highly excited initial states {v = 10,40,62 out of 66 bound levels). For moderately excited 
vibrational levels, a direct pathway from the initial to the target state exists, and a single central frequency in the 
guess field was sufficient to obtain an optimal field. For higher excited levels, a multistep scheme is required since 
initial and target state have good Franck-Condon overlap with different vibrational levels of the excited state. This 
was accounted for by assuming two central frequencies in the guess field. Experimentally, pulses with two central 
frequencies can be realized by sending a pulse with a single central frequency through a parametric amplifier. 



14 



0.08 



0.04 - 







0.08 



0.04 



- 1 1 1 1 

fr.\ J 
W , 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 - 




v=61 




v=63 - 


.. f 1 ■ 

i ' 

m ' ' 

Mil 

-J V: 




/\: 

^1 1*1 


1. 



T — I — I — I — I — I — I — I — r 
(b) 



"1 — I — I — I — rrn — I — I — r 



v=22 

V=23 

V^21 





0.04 







1 — rn — I — I — I — I — r~r~i — i — i — r~r 



V =105 

V =104 

V=103 

V =106 



o.oe 



T — rn — I — r~m — m — i — i — r~i — i — i — i — i — r~r 

hCd) 





r ' J ■ I I I >-'-r- 



- 0.08 



0.04 







— I — I — I — I — I — I — m — r~i — I — I — I — r"T ■( '"i" i— j 



V"=107 
V"=108 



V"=102 H 

v =ioi 



T — ( — I — I — I — I — I — I — I — I — I — I — I — I — I — I — ] — I — T 



V=51 

V'=9 

V"=8 

V =7 



- 0.2 



I I r- 



0.1 







§ 0.5 1 1.5 

time [ps] 



20 



0.5 1 1.5 

time [ps] 



FIG. 8: Vibrational level population of the dynamics with the optimal field for the electronically excited state (c)-(f) and the 
ground state (a)-(b) 



The two control knobs of the optimization algorithm are the pulse intensity and the optimization time. A compar- 
atively high intensity is needed to excite several quantum pathways whose interference constitutes the solution. The 
optimization time is usually chosen to correspond to the longest vibrational period of the system. We have explored 
ways to restrict both intensity and optimization time. This was possible by using information from optimal fields 
of a previous run of the algorithm with large intensity / long optimization time in the guess field which was used 
as input in a second run of the algorithm with smaller intensity and shorter optimization time. Optimal fields with 
integrated pulse energy on the order of milliJoule and overall duration on the order of a few picoseconds have thus 
been obtained. 

The optimal control task we have investigated in this paper is much more difficult than those of previous control 
studies on the sodium dimer, such as populating certain ionization channels jsJl or displacing the ground state 
Gaussian wave packet [3^j. The fact that solutions could still be found, is due to the modification of the optimal 
control algorithm to calculate the change in field rather than the new field itself |33| . This ensures a smooth convergence 
of the field toward the optimal one and renders the algorithm very powerful. 

We have neglected several possible loss channels. Two-photon transitions could populate an autoionizing state. 
However, since molecular autoionization occurs at small distances, the cross-sections are very small l4^. The field 
could moreover transfer angular momentum to the system thus exciting rotations. Both loss channels can be modelled 
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by adding one or a few electronic states to the model. In general, additional channels do not necessarily pose a 
problem, they may even lead to more participating pathways improving the control jsj]. Intuitively, the field is 
able to avoid population of these loss channels if its phase is orthogonal to the transition dipole. It has been shown 
in local optimization calculations |U |4^, that this condition on the phase of the field can be used to lock the 
population in the desired channels. In global optimization, convergence can be facilitated by imposing a penalty on 
the loss channels Multiphoton ionization of the sodium atoms is more dangerous than excitation of rotations 

and molecular autoionization. It is rather likely to occur at the intensities we have found. Our calculations should 
therefore be repeated by adding an A^-level Hamiltonian describing the atomic levels to the model. The task of the 
modified control problem is to transfer the vibrationally excited molecules to v = while avoiding population of the 
ionization loss channels at any time during the optimization interval. The latter has to be formulated as an additional 
constraint. The function g(e) does then also depend on the state (p, g{e, ip), and the equation for the new field has to 
be rederived. Work in this direction is in progress. 

In a recent experiment j5| , ultracold Na2 molecules have been produced from an atomic sodium condensate with the 
molecules occupying the i; = 14 vibrational level of the a^I]+ state. This vibrational level is very weakly bound (< 1 
cm~^), and at such small binding energies the hyperfine interaction mixes singlet and triplet states Note that 
in the asymptotic region the wave function of the triplet v = 14 level is identical to that of the singlet w = 64 level. 
The coupling between hyperfine states implies that also vibrational levels of the singlet state are populated, hence 
our results considering the X^E+ should be applicable to the MIT experiment. Second, to obtain quantitative results 
our model should be improved to include all the states which are coupled by hyperfine interaction. This requires, 
however, several technical improvements and is therefore the subject of future work. 

The MIT experiment raises another question regarding the feasibility of optimal control experiments with ultracold 
molecules. The reported lifetime of the molecular cloud is a few ms 0- Conventional pulse shapers operate at about 1 
kHz resulting in one cycle of the learning loop per millisecond. In standard control experiments a few thousand cycles 
can therefore be performed in a comparatively short time. If the MIT experiment were to be combined with a control 
scheme, only one or two cycles could be performed before the molecules are lost from the trap. It therefore seems 
necessary to provide more information from theoretical considerations than simply the required spectral bandwidth 
and intensity of the pulse. A promising approach has been suggested by de Vivie-Riedle and coworkers 51] who 
defined the mask pattern of the optimal pulse by discrete Fourier transforms. Provided the calculated pulse is not 
too complex, this mask pattern can be directly fed into the pulse shaper avoiding the many cycles of a learning loop. 

In addition to accounting for ionization losses and hyperfine interaction, our two-state model should be improved to 
consider rotations. This is particularly important because the molecules are not necessarily created with J = in the 
experiments using Feshbach resonances or three-body recombination. Since the rotational excitation is rather small, 
it will be sufficient to treat it by additional electronic states, for which the centrifugal barrier for each J is added to 
the atom-atom interaction. One may speculate that some rotational excitation may even facilitate the stabilization 
process since the initial wave function due to the centrifugal barrier is more localized than the one for J — 0. 

Another possibility which leads to a more confined initial state for the stabilization process is given by forming the 
molecules by photoassociation. Two-color photoassociation or spontaneous emission after one-color photoassociation 
can directly populate levels around u = 40. We have shown in our calculations that the optimal control task for such 
vibrational states which are still highly excited but already well-bound is much easier than for the last bound levels. 
The required intensities are hence more moderate (Cf. an integrated pulse energy of 1.5 mJ for v = AO vs. 4 mJ for 
V = 62). 

Finally, our calculations prompt the question whether spectrally simpler fields than those obtained in this study can 
be guessed. Our population analysis has shown that at least two cycles through the excited state are required. The 
first cycle populates excited state vibrational levels which have good Franck-Condon overlap with the initial state, 
and the second populates excited state vibrational levels with good Franck-Condon overlap with the target state. One 
could hence divide these steps, optimizing each separately. First calculations show that the intensities required for one 
such step are much lower than the ones completing the control task at once. Note that while the intensity is reduced, 
the overall time (including all steps) is prolonged. Another possibility is given by finding a solution for each of the 
steps intuitively without the optimal control algorithm, for example chirping the sequence of pulses. Furthermore, 
stabilization via optimal control could be combined with chirped pulse photoassociation [s^ . l53| . For the Cesium 
dimer, it has been shown that a spatially focussed wave packet can be formed on the excited state [s^. To stabilize 
such a state to the ground vibrational level of the electronic ground state via OCT should also require comparatively 
little intensity. Our results furthermore indicate that STIRAP or any other adiabtic process is not likely to be a 
successful tool in transferring the highly excited molecules to the vibrational ground state. The multistep scheme 
implies that there is no direct adiabatic path, i.e. the Franck-Condon overlaps of a direct path are extremely small. 
Correspondingly, the intensity required for a direct adiabatic path is expected to be orders of magnitude larger than 
that obtained by us for an optimal field. Adiabaticity moreover requires the transfer to be slow. If we quantify slow 
as being at least one order of magnitude slower than the longest vibrational period of the problem, the process has to 
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occur on the time scale of nanoseconds in the case of Na2 and even longer times for heavier molecules. These are the 
time scales where loss processes become important and where condensate dynamics and molecular dynamics become 
entangled. 

In conclusion, we have applied in this study optimal control theory to the stabilization of ultracold molecules. 
While we have started with a comparatively simple model, our present results point toward several directions of 
improvements which in turn correspond to different experimental schemes. We thus hope that the present work 
stimulates further research uniting the fields of cold molecules and optimal control. 
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APPENDIX A: BRIEF REVIEW OF OPTIMAL CONTROL THEORY WITH THE KROTOV METHOD 

The Krotov method is one of several methods to derive an iterative algorithm which connects the objective of 
optimal control, calculated at the final time t = T, with the knowledge of the state of the system at intermediate 
times < i < T. A general review was given in Ref. |i22il] . while a review for linear problems in the density matrix 
formalism was reported in Ref. ,37,J . In this paper, we need an even simpler version, that is the Krotov method for 
linear optimization of state-to-state transitions which we sketch in the following. 

The equation of motion of the system is the Schrodinger equation, 

^|(p(t))^-lH(e)|^(t)), 

whose right-hand-side is abbreviated by 

|/(t,<^(i),e)) = --^H(£)|^W). (Al) 

We omit the ket notation when just indicating the dependence on \ip{t)). The starting point is now the modified 
objective J of Eq. jSJ, which is a functional of the state of the system and the field, J[(p, e]. Due to integral term in 
Eq. (O, J is a functional of the system's state at all times, i.e. a functional of the system "trajectory". The problem 
is now to find conditions on J which maximize the objective F. 
To this end, one introduces an auxiliary functional L[ip, e; $], 

L[^,e;<i>] =G(^(r))-<i>(0,¥P,)- / Rit,vit),e)dt (A2) 

Jo 

with 

G(^(T)) = -F(^(T)) + ci>(T, ^(T)) , (A3) 



\ dip 



(A4) 



and $ = ^{t, If) an arbitrary continuously differentiable real function j23,l23- Recall that jtpi) = |</j(< = 0)), and \(pf) 
is the target state. Note that the derivative of $ with respect to the state means derivative with respect to both real 
and imaginary part of \-^) — | + ^^fj)- If the state \ip(t)) and the field e(t) obey the Schrodinger equation, 
R is given by i? = + It can then be shown 22] that L[ip, e; $] = J[ip, e] for any scalar function i.e. the 

minimization of J is equivalent to the minimization of L. In the latter one has complete freedom in the choice of $. 
This property is used to construct an iterative algorithm. One first chooses $ such that L is maximum with respect 
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to \ip{t)), i.e. the worst case. A new field is then derived from the condition of maximizing R which in turn leads to 
minimization of L. 

The first step (maximization of L) is equivalent to maximizing G while minimizing R. This step is taken at the old 
state, |</?(t))". To simplify the equations, the maximum and mininum conditions are relaxed to extremum conditions. 
An additional condition in the end has then to ensure that the algorithm indeed improves the objective in each 
iteration. The extremum conditions are written as 
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dip 
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The partial derivatives are again with respect to both real and imaginary part of Eq. ljA5|l leads to 
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Using that |/) = d/ dit\ip) and calculating \d/d'p){f \ from Eq. ljAl|l . we obtain 
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Eq. ljA6|l together with the definition of G{(pf) (Eq. 
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(A5) 
(A6) 



(A7) 



(A8) 



Recall that F is the original objective, and \dF/dip) can thus be explicitly calculated. Eq. ljA7|l can then be interpreted 
as first order differential equation for \j{t)) — f^(i)/ with the "initial" condition (at time T) given by Eq. jSHl- 
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The knowledge of Eqs. (|A7|) and (jA8|) is therefore sufficient to determine ^{t,ip) to first order in \ip{t)), 
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From $(1), G and R can be constructed to first order, G^^^ and i?'^^^. Thereby the first step of the algorithm for a 
linear optimization problem is completed. Note that <i>^^^ is never explicitly calculated, the equation of motion for 
\^it)) (Eq. HA7|l ) will be solved instead. 

To realize the second step, R^^^ must be maximized with respect to the field. This step is supposed to lead to a new 
field, e""*"^, which generates a new "trajectory", |(y3(i)}"+^, according to the Schrodinger equation. We again relax the 
maximum to an extremum condition, 



which gives 
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and therefore 



dg{e) 



de 



de 



(All) 



With our choice of the constraint g{e), Eq. and using H = Hq + fieit) in the calculation of \df /de) , we obtain 
for the new field 



(A12) 



The time evolution of the new state is given = U(i, 0; £"+^)|(/5i}, and \j(t)) is obtained by formally solving 

Eq. jSU, |7(t)) = 0(i,r;e")|7(T)}, 



e"+i(t) =£"(t) + 



2a 



3m 



(7(T)|U'(i,T;£")AU(i,0;e"+0|(p.) 



Using the definition of the objective F, Eq. Q), |7(T)) is determined according to Eq. (|A8I) . |7(T)) = c|(^/) with the 
coefficient 



c=(^/|^(r)r = (^/|U(T,0;e")|^.) 



(A13) 



and finally Eq. {Tj) is obtained. 

Since we relaxed the maximum and mininum conditions to extremum conditions, we need to ensure that each 
iteration indeed improves the objective, J[(/3",e"] — J[(^"+^; £"+^] > 0, 



Ai+ /" A2{t)dt 
Jo 



with 



A sufficient condition for J[y" ,£"] - J[(^"+\£"+i] > is Ai, A2(t) > 0. Constructing G^^^ from Eq. ijUl), and using 
|7(r)) = c\ipf) with Eq. PT3|) we see that 



Ai 



> 0. 



(AM) 



(<^,|U"(r,0;£")-U"(T,0;£"+i)|^/) 
Because of the linearity of f{t, tp, e) w.r.t. \ip), i?(^^(t, ip, e") = — (?(£") for any \ip), and we find 

A2 (t) = [i, </.(<)"+! , £"+!] - [t, ip{tr+' , £"] . 

Constructing i?*^^^ from Eq. l|A9p and using again the linearity of the equation of motion, we obtain as condition for 
monotonous convergence 



A2(i) = -.g(£^)+ff(£°) + (£!-£") If 



> 0. 



(A15) 



Eq. ljAll|l and Eq. ljA15|l are the central equations of the algorithm. The constraints which we formulate in g{e) 
must fulfiU Eq. ![KT^ to ensure convergence of the algorithm. Throughout this work wc have used the constraint of 
restricted change in pulse energy, i.e. g{e) as given by Eq. lO, which respects Eq. IjAlSj) for all a as one can easily 
verify. 



APPENDIX B: RESTRICTION OF THE SPECTRAL BANDWIDTH 



We tried other functional forms of the constraint g{e), in particular 



gie) = aiit) - e"y - a2{t) {e ~ e'^') 



(Bl) 
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with ai = aoi/S{t), where the second term restricts the new field to a reference field e''®^ with a prespecified bandwidth 
(note, that the negative sign of this term implies maximization, i.e. the new field should be as close as possible to the 
reference field). Condition l)A15|l is fulfilled for all aoi > q;o2 > 0. We found no effect for q;o2 ^ ctoi as not enough 
weight is put on the second constraint. Still 96% of the objective can be reached for ao2 close to aQi, but the spectral 
bandwidth is slightly reduced. In particular, spurious high-frequency components in the optimal puis can be avoided. 
The overall restriction of the bandwidth, however, appears to be insufficient. 

In fact, it is not surprising that a condition formulated in time domain as Eq. (|B1|) cannot achieve the desired 
control over a frequency domain property. This is a general problem of global (in time) optimal control schemes. 
Attempts to restrict the spectral bandwidth of the puis have been made before However, this is possible only 

at the cost of monotonic convergence. An approach which unifies global optimal control with constraints in frequency 
domain will have to treat time and frequency on the same footing. This shall be the subject of future work. 
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